#include<list>
#include<vector>
#include<boost/shared_ptr.hpp>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#ifdef USE_MPI
	#include<mpi.h>
#endif

#define EPS 5.96e-08 //floating point machine-epsilon

class cListElement;
typedef unsigned int uint;
typedef boost::shared_ptr<cListElement> cListElementPtr;
class cLinkedList;

class cOpenCLInterface;

class cFlySim{
public:
  static void initialize(const char * xmlfile);
  static int get_rank();
  static int get_ranks();
  cFlySim();
  ~cFlySim();
  static string analysis;
protected:
  // The following are essentially global variables  
  static bool is_initialized;
  static bool enable_gpu;
  // user specific variables
  static const float SQRT_2;
  static long int ranseed;
  static string outputfile;
  static int gciTotalNumberOfFlies;
  static int gciTotalNumberOfMales;
  static int gciNumberOfPetriDishes;
  static int debugSerialMPI;
  static float dThreshold;
  static float dMaleAveAve;
  static float dMaleAveSD;
  static float dMaleVarAve;
  static float dMaleVarSD;
  static float dFemaleAveAve;
  static float dFemaleAveSD;
  static float dFemaleVarAve;
  static float dFemaleVarSD;
  static float dMFCovarAve;
  static float dMFCovarSD;
  static float dMFCovCorrection;
  //the following are the 'new' stats
  static float dMpMAve;
  static float dMpMSD;
  static float dFpMAve;
  static float dFpMSD;
  static float dMpFAve;
  static float dMpFSD;
  static float dFpFAve;
  static float dFpFSD;
  //the following are the new weighting coefficients
  static float dAvBoysWeight;
  static float dAvGirlsWeight;
  static float dVarBoysWeight;
  static float dVarGirlsWeight;
  static float dCovBoysGirlsWeight;	  
  static float dMalesPerMaleWeight;
  static float dFemalesPerMaleWeight;
  static float dMalesPerFemaleWeight;
  static float dFemalesPerFemaleWeight;
    
  //the number of stats you want to store
  static int ciNoStats;
  static float dSamplingInterval;
  static int iNSamples;
  static int iSamplingInterval;
  static int iTimeToStartSampling;
  // ABC parameters
  static bool bDoABC;
  static int gciNoOfIterationsInABCAnalysis;
  static int gciBurnInPeriod;
  static int gciABCSamplingPeriod;
  static int gciNoOfSimulationsPerProposal;
  static float dABCTolerance;
  static float dABCNewTolerance;
  static float dNewSumsSquaresValue;
  static float dSquaredDistance;
  static bool bUseBortotScheme;
  static bool bToleranceChanged;
  static float cdProbChangeTolerance;
  static float cdProbToleranceMovesDown;
  static float gdOldBortotTolerance; // book-keeps the value of the threshold when we are changing it
  static float gdMaxTolerance;
  static float gdMinTolerance;
  static float gdMaxToleranceChange;
  static float gdInitialToleranceWhenUsingBortot;
  static float gdPriorLambdaForTolerance;
  static int gciNumberOfFlyParameters;
  static bool bOutputSummaryStats;


  // OpenCL interface
  static cOpenCLInterface * clInterface;
  // The GSL random number generator
  static gsl_rng * rg ;
  // The outputfile.  Any class can write debugging output
  static ofstream ofile;
#ifdef USE_MPI
  static int mpi_rank, mpi_ranks;
#endif

};

// We need to tell the next class (cPetriDish) that there is another class that exists and that it is called 'cFly')

class cFly; 
  
class cPetriDish:public cFlySim  // This is where flies hang out 
{
public:
  cPetriDish(); // the static constructor - it sets up new Petri dishes
  ~cPetriDish(); 
  friend class cExperiment;
  friend class cFly;
private:
  int iIndex; // A unique identifier for this petri dish
  float dFoodQuality;  /// how good is the food here?
  int iHowManyMalesAreHere; // To save time, we also keep track of how many of the flies here are male
  int iHowManyFemalesAreHere; // To save time, we also keep track of how many of the flies here are female

  //bool bCheckOurRecordsAreGood(); // This function just checks that we have the correct number of males and females recorded as being here
  void PrintDetailsOfThisDish(); // Prints out details of this dish
};

class cache_container_t{
public:
  cache_container_t(int iTotalOthers);
  ~cache_container_t();
  float * dAttractionConstant; // array of 4 interaction types
  float * dPreComputedPowers; // matrix of 4 inter types * number of partners
  float * dPreComputedInteractions; // interactions between M to F precomputed powers
};

class cFlyCategory:public cFlySim{
public:
  cFlyCategory();
  ~cFlyCategory();
  static const int INTERACTION_MM = 0;
  static const int INTERACTION_MF = 1;
  static const int INTERACTION_FM = 2;
  static const int INTERACTION_FF = 3;
  void commitParameters();
  void rollbackParameters();
  void setMaleAttractionToMales(float); 
  void setMaleAttractionToFemales(float);
  void setFemaleAttractionToMales(float); 
  void setFemaleAttractionToFemales(float); 
  friend class cFly;
private:
  // will get called by updateAttractionConstant if necessary
  void updatePreComputedPowers(int interactionType);
  // Cached computed powers
  cache_container_t * current_cache;
  cache_container_t * previous_cache;
  //float * dPreComputedPowers; // matrix of 4 inter types * number of partners
  //float * dPreviousPreComputedPowers; // matrix of 4 inter types * number of partners
  //float * dPreviousAttractionConstant; // array of 4 interaction types
  //float * dAttractionConstant; // array of 4 interaction types
  float * dPowers;
  int iTotalOthers;
  void copyCache(cache_container_t * in, cache_container_t * out);
  //float MaleAttractionToMales; 
  //float MaleAttractionToFemales; 
  //float FemaleAttractionToMales; 
  //float FemaleAttractionToFemales; 
};

class cFly:public cFlySim // This is a fly
{
public:
  cFly(cFlyCategory * flyCategory); // static constructor for new flies
  ~cFly(); 
  friend class cExperiment;
  float getMaleAttractionToMales(); 
  float getMaleAttractionToFemales();
  float getFemaleAttractionToMales(); 
  float getFemaleAttractionToFemales(); 
private:
  int iTotalOthers;
  cFlyCategory * flyCategory;
  char cSex; // What sex am I?
  cPetriDish *WhereAmI; // A pointer to the petri dish that the fly currently calls home 
  void Print(bool); // prints detail of this fly; if bool==1 do it on a new line
  float dLeaveRate; // The rate at which I leave my current location given who is currently around. It will depend on the above 8 parameters.
  int LandOnAPetriDish(); // A fly decides to land on a petri dish
  int LeaveThisPetriDish(int); // A fly decides to leave a petri dish
  int iID; // a unique identifier for this fly
  int iGenotype; // what genotype am I?

  //Fly parameters
  
  float RateLeaveMale;
  float RateJoinMale;
  float RateLeaveFemale;
  float RateJoinFemale;
  float expRateLeaveMale;
  float expRateJoinMale;
  float expRateLeaveFemale;
  float expRateJoinFemale;

  void dOldCalcRateOfLeaving(); // This function calculates the rate at which a fly wishes to pack its bags and leave its perti dish
  
  void dCalcRateOfLeaving(); // This function calculates the rate at which a fly wishes to pack its bags and leave its perti dish
  float dRateByInteraction(int interactionType);
  // More efficient way of calculating rates that we often recalculate
  //float dRateMM(float dAttraction, int iHowManyMales);
  //float dRateMF(float dAttraction, int iHowManyFemales);
  //float dRateFM(float dAttraction, int iHowManyMales);
  //float dRateFF(float dAttraction, int iHowManyFemales);

};

//********* list classes

class cListElement:public cFlySim // our linked list is going to consist of a list of these things all linked together by pointers
{
public:
  cListElement(); // static constructor
  cListElement(const cListElement & listElement); // static constructor
  ~cListElement(); 
  friend class cExperiment;
  friend class cLinkedList;
private:
  float * dStats; // an array to store the statistics
  cListElement* pNext; // points to the next element in the linked list
  cListElement* pPrev; // points to the previous element in the linked list

  // functions
  void Print(); // prints this element
  void ListWrite();
};


class cLinkedList:public cFlySim// this is the linked list
{
public:
  cLinkedList(); // static constructor - this sets the thing up the first time you create it
  friend class cExperiment;
private:
  void WriteTheEntireList(); // prints each element of the linked list
  void EmptyTheList(); // empties and deletes the list when you are done with it (otherwise it may live in memory for a long while)
  void AddAnElement(); // adds an element to the end of the list
  void RemoveAnElement(); // Remove the last list element
  cListElement BottomOfList; // This is the 'anchor' for our collection of list elements
  cListElement* pTopOfList; // this points to the last element in the list
  
  // We need a couple of variables to help make our life easier
  cListElement BottomOfPileOfSpareListElements; // It pays to keep list elements you are not 
  // currently using rather than deleting them and then having to redeclare them when you want more (that takes a while)
  
  cListElement *pTopOfPileOfSpareListElements;

  // functions    
  cListElement* GetAnElement(); //  gets a new list element
  void PrintTheEntireList(); // prints each element of the linked list
};

const int ciNoOfParameters=6;
const int ciDimensionOfData=ciNoOfParameters;
const int ciNoOfAcceptancesPerGeneration=500;//1000;


class SMC:public cFlySim
{
public:

    double TestParamValues[ciNoOfParameters];
    double StdDev;
    const static int ciNoOfTestDatasets=10;
    double dRejectionThreshold;
    int iMaxNumberOfAttempts;
    int iNumberOfGenerationsForSMCAnalysis;
    int ciNoOfAcceptsNeeded;
    double TestData[ciNoOfTestDatasets][ciDimensionOfData];

    double SampledParameters[ciNoOfParameters];
    double GeneratedData[ciDimensionOfData];
    double NewAcceptedParameterValues[ciNoOfAcceptancesPerGeneration][ciNoOfParameters];
    double OldAcceptedParameterValues[ciNoOfAcceptancesPerGeneration][ciNoOfParameters];
    double dWeightsOfOldParametersValues[ciNoOfAcceptancesPerGeneration];
    double dWeightsOfNewParametersValues[ciNoOfAcceptancesPerGeneration];
    float dParameterLowerBounds[ciNoOfParameters];
    float dParameterUpperBounds[ciNoOfParameters];


  SMC();
  ~SMC();
  void InitializeArrays();
  void PrintParameterValues(bool);
  void PrintWeightedParameterValues();
  void RunSMCAnalysis();

//  static const int ciNoOfAcceptancesPerGeneration=1000;
	

private:
  //ofstream ofile("SMCout1.txt",ios::out);
 // static const int ciNoOfParameters=3;
  //float dParameterLowerBounds[ciNoOfParameters];
 // float dParameterUpperBounds[ciNoOfParameters];
  //static const int ciDimensionOfData=3;
  //ustatic const int ciNoOfAcceptancesPerGeneration=1000;
//  static const int ciNoOfAcceptancesPerGeneration;
 // int iNumberOfGenerationsForSMCAnalysis;
  float ThisGensParams[ciNoOfAcceptancesPerGeneration][ciNoOfParameters];
  int TempArray[ciNoOfParameters];
  static const bool bDoPureRejection=0; // generally this will be set to 0, but if we set it to one it will do rejection without SMC ( y setting the number of generations equal to 1

  //float TestParamValues[ciNoOfParameters];
  //static const float StdDev=0.5;
  //static const int ciNoOfTestDatasets=10;
  //static const float dRejectionThreshold=12;
  //static const int iMaxNumberOfAttempts=1600000;
 // static const int ciNoOfAcceptsNeeded=ciNoOfAcceptancesPerGeneration;
  /*float TestData[ciNoOfTestDatasets][ciDimensionOfData];
  float SampledParameters[ciNoOfParameters];
  float GeneratedData[ciDimensionOfData];
  float NewAcceptedParameterValues[ciNoOfAcceptancesPerGeneration][ciNoOfParameters];
  float OldAcceptedParameterValues[ciNoOfAcceptancesPerGeneration][ciNoOfParameters];
  float dWeightsOfOldParametersValues[ciNoOfAcceptancesPerGeneration];
  float dWeightsOfNewParametersValues[ciNoOfAcceptancesPerGeneration];
*/

  void BuildTestData();
  void PrintTestData();
  void DoThisGenerationsRejection(int);
  void GenerateParameters(int);
  void GenerateData();
  float MeasureDistanceBetweenDataAndParameters();
  void UpdateParameterValuesAndWeights();
  void NormalizeNewParameterWeights();
  void CalculateParameterWeights(int);
  float SampleFromParameterPrior(int);
};


class cExperiment:public cFlySim  // A collection of petri dishes and flies
{
public:
  cExperiment(); // static constructor
  ~cExperiment();
  void run();
  
private:
  vector<cFlyCategory *> flyCategories;
  vector<cFly> Flies; // these are all the flies
  vector<cListElement> statsVector;
  float * dStatsVec;
  int iTotalTimeOfExperiment;
  float p;
  cPetriDish * Dish; // these are all the petri dishes - by convention we will make petri dish #0 represent not being on a petri dish
  //float PetriDishSummaryStats[gciNumberOfPetriDishes+1][1]; // stores summary statistics of interest about the petri dishes
  float * PetriDishSummaryStats;
  //float TargetPetriDishSummaryStats[gciNumberOfPetriDishes+1][1]; // This stores the target (i.e. observed) values of the statistics
  float * TargetPetriDishSummaryStats;
  long int iSamplesTaken; // records how many times we sample from the population
  float dTime; // records for how long the experiment has been running
  float dLastSampleTime; // records the time at which we last sampled the summary stats
  long int iMoveCounter; // counts the total number of occasions at which any fly has decided to move - just out of interest
  long int iNumberOfAcceptances; // Counts how many ABC moves have been accepted
  float * dOldParameterValues; // stores the old parameter values for fly movement so that we can return to them in an ABC move is reject
  int iTimesAFlyMoved;
  int iTimesAMaleMoved;
  int iTimesAFemaleMoved;
  // PM added these delcarations:
  int giPetriDishLeft,giPetriDishArrivedAt;
  // function prototype for reading configuration settings from XML file
  void OldRandomSetup(); // Drops flies randomly onto petri dishes
  void RandomSetup(); // Drops flies randomly onto petri dishes
  void PrintCurrent(); // outputs current details of the experiment
  void OldNextEvent(); // works out what the next event is and implements it 
  void NextEvent(); // works out what the next event is and implements it 
  void PackUpAndLeaveTown(int); // The fly indexed by the argument leaves its current location and moves somewhere else
  void UpdatePetriDishes(cFly*,int,int); // updates the petri dishes that the fly is moving from and to
  void PrintSummaryOfPetriDishes(); // output how many flies there are on each petri dish, including the 'outside environment', which is listed first)
  void UpdateSummaryStats(); // updates the summary stats we are keeping
  void UpdateList(float); // new function to update the values in the list
  void OutputSummaryStats(); // outputs the summary stats - NOW updates and outputs
    void OutputSummaryStatsAtEndOfSim(); // outputs summary stats when we are not doing ABC - PM added on 1/15/13
  void CheckState(); // this function just checks that all our recrod-keeping variables are up-to-date. It's for debugging purposes
  void Evolve(int,bool,bool); // Evolves the population through a period of time of length 'float'; if 1st 'bool'==1 we record and output summary stats; the 2nd 'bool' controls whether those summary stats are output to screen
 float CalcEventRate(); // Calculates the rate at which the next event occurs

  // the following functions are to do with the ABC analysis
  void ResetAllPetriDishStatistics(); // initializes everything to 0, ready for the next simulation
  void ABC(); // performs a time-saving ABC analysis
  void oldABC(); // performs an ABC analysis
  void InitializeEvolutionaryParameters(bool); // chooses a set of parameter values to start from
  void OutputParameterValues(int,bool); // outputs the parameters of the flies
  bool IsSimulatedDataWithinTolerance(); // checks whether each summary stat from the simluated data is within 'dABCTolerance' of the observed value
  void ProposeNewState(); // proposes a new set of parameter values
  void ReturnToOldState(); //discards the new set of parameters and returns to the old ones if the move is rejected
  void ChangeTolerance(); // PM added this to change the threshold automatically  
  float CalculateHastingsRatioForToleranceChange(); // PM added this to calculate the HAstings Ratio for proposal changes
  bool BortotIsSimulatedDataWithinTolerance();  // thius function works differently when we are using Bortot scheme
  float PriorDensityForTolerance(float); // This function defines the prior for the tolerance
};

